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ABSTRACT 

The point charge model is used to calculate the crystal structure of sulfuric 
acid (H2SQ04) with the 6-31G** basis set. The point charge model accurately 
reproduces the structural trends which occur in transforming from the gas to the 
solid phase. The calculated crystal structure of sulfuric acid is in reasonably good 
agreement with both the X-ray and neutron-diffraction structures. 

The point charge model is shown to precisely simulate the deformation 
forces which are present in the solid upon crystallization. The point charge model 
exhibits a definite shift of electron density from the bridging hydrogens to the 
acceptor atoms, identical to those found in other ab initio studies. The calculated 
crystal structures are insensitive to the magnitude of the point charges. 

Two new iterative techniques, using the point charge model, are introduced 
which give superior results to any of the single optimization cycle methods. These 
iterative techniques account for any forces which are not electrostatic by nature. 

The gas phase structure of H2SQOgq is optimized using the STO-3G*, 4- 
31G**, and 6-31G** basis sets. Comparisons are made between the programs 


TEXAS, TEXAS(D), and Gaussian-82. 
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CHAPTER ONE 
INTRODUCTION 

Over the years, many researchers have shown the ability to apply quantum 
theory to chemistry in order to supplement experimental results. However, it has 
only been in the past decade with the great advances in computers coupled with the 
development of efficient algorithms that ab initio calculations have become practical 
on molecules of more than a few atoms. 

The term "ab initio" is frequently cited in the literature but rarely defined. 
Ab initio calculations can basically be thought of as those which treat all the 
electrons and do not utilize parameters which are fit to experimental data as 
compared to semi-empirical calculations. For example, a Hartree-Fock calculation 
of a restricted closed-shell single determinant wave function is one in which no 
approximation is made to the integrals or the electronic Hamiltonian. However, it 
is completely specified by the choice of a basis set and the coordinates of the 
nuclei.! 

One area in which ab initio calculations have been particularly successful is 
in evaluation of molecular structure. The theoretical determination of a molecular 
structure involves minimizing the total energy of the molecule with respect to 
simultaneous variation of the internal coordinates.2 The application of this 
procedure to complex molecules became practical only after the introduction of the 
gradient method by Pulay. This method greatly reduced the time to calculate an 
optimum geometry especially for large molecules.3 The resulting optimized 
structure has a known basis set dependance which results in a constant error for a 





given parameter called the offset value. This offset value is reasonably constant 
over a wide range of molecules.4 Boggs and coworkers have reported extremely 
accurate structures by combining experimental and theoretical data. The structures 
presented in their study are more accurate than either theoretical or experimental 
techniques could produce alone.? Not only have ab initio calculations been able to 
accurately reproduce experimental bond lengths, bond angles, and conformations 
but in some cases ab initio calculations have actually shown where experiments 
have determined incorrect structural parameters. Although the latter is the 
exception, it shows both the value and validity of ab initio calculations as a 
supplement to experimental data. 

The ability to determine accurate molecular structures is the key to many 
areas of chemistry. The structure of a molecule not only defines its physical 
properties but also influences how it reacts with other molecules. This information 
is vital to such diverse areas as the development of "stable", powerful explosives 
and propellants, the synthesis of effective catalysts, the understanding of how 
enzymes with only specific structures can fit into certain substrates, the 
determination of complex reaction mechanisms, and the production of 
superconductors. However, to a chemist, the most fundamental concern which 
structural information can shed light on is chemical bonding.© The understanding 
of chemical bonding is fundamental to all areas of science. 

One area in which there is a great deal of interest is the structural differences 
which occur in some molecules. These differences can be between the same type 
of bond (eg, C-H) or functional group (eg, CH3) over a variety of similar 


compounds with different substituents. These changes can also occur within the 
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same molecule when it is in a different phase or chemical environment (eg, in 
solution). These structural differences need not be large to be important. In fact, 
they can be very small. Differences on the order of 0.01 A in a bond length or a 
couple of degrees in a bond angle may reveal extremely valuable structural 
information if they are reliably determined. However, only real structural 
differences should be considered and any apparent differences which are detected 
must be eliminated if meaningful comparisons are to be made.’ Apparent 
differences can arise from either the particular interaction of radiation with matter or 
because of the vibrational averaging which occurs in a particular experiment. It 1s 
important to note that all experimental techniques produce geometrical parameters 
which are based on averaging of the molecular vibrations. However, at least 
theoretically, ab initio calculations give equilibrium intermolecular distances which 
correspond to the minimum of the potential energy surface and are not subject to 
molecular vibrations. 

The vast majority of all ab initio calculations are performed on isolated 
"gas" phase molecules. The reason for this is quite simple - computer capacity. 
The most accurate results come from ab initio calculations performed with large 
basis sets. However, current computer capacity limits large basis set calculations to 
no more than about two dozen atoms. This limit of two dozen atoms with a fairly 
large basis set is pushing the capacity of even the Cray supercomputer. Thus, in 
order to conduct conventional ab initio calculations on solids, the size of the basis 
set would have to be greatly reduced, a lower level of theory must be used, or 
assumptions must be made about the symmetry of the molecule. These restrictions 


almost always result in less accurate molecular parameters. 
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The intent of my work was to study the structural changes which occur in a 
molecule between the gas and solid phase through the use of ab initio calculations. 
In selecting a particular molecule to study, I established several criteria. First, the 
number of atoms in the molecule must be small enough that a large basis set could 
be used to give accurate results but be within the capacity of the Cray 
supercomputer. All of the calculations in this study were performed on the Cray X- 
MP/24 at the University of Texas Center for High Performance Computing 
(UTCHPC). Second, there must be an accurate experimental determination of the 
gas and solid phase structure so that comparisons could be made with the 
calculations. Finally, if possible the molecule should be one which is of general 
interest. After careful consideration, sulfuric acid (H2SO4) was selected as the 
molecule to study which met all of the above criteria. 

In order to reduce the calculation of solids to one that can be handled by 
present computers, a simpler method must be used in order to simulate the effects 
of the crystal lattice. The problem is basically how to accurately represent the 
molecular-charge distribution of the surrounding molecules in the solid. The 
simplest is to use the Mulliken population analysis.8 The Mulliken population 
analysis comes from parameters that are already calculated in quantum mechanical 
calculations. The major problem with this method is that the assignment of charges 
to the atoms is rather arbitrary and is often very basis set dependent.8 However, a 
previous study done using the Mulliken population analysis and point charges was 
able to accurately reproduce the crystal structure of cyanoformamide 
(NCCONH)2).? The purpose of my study was to apply the point charge method to 


a different type of molecule and to develop some new techniques in order to see 
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their effect. The theoretical methods used in this study are presented in Chapter 2. 
Chapter 2 includes not only the basis set and level of theory considerations but also 
a review of the various methods which have been used to simulate crystal field 
effects. Chapter 3 contains the sulfuric acid results including the optimized gas 
phase geometry and the various techniques used to calculate the solid geometry. 


_ Chapter 4 presents my conclusions and suggestions for future studies. 





CHAPTER TWO 
METHODS 

In any ab initio calculation, the selection of the basis set and the level of 
theory are crucial. A decision must be made between the computational cost and 
the desired accuracy. Obviously, the most accurate results are desired from any 
calculation. However, a structure accurate to 0.01 A at one-hundredth the cost of 
one which is accurate to 0.001 A maybe preferred if the first structure demonstrates 
the effects which are being studied. Although, this is a hypothetical example, it 
points out the importance of extremely careful selection of the basis set and level of 
theory. In my particular study, the selection of a method to simulate the crystal- 
field effects is also a major factor. Again, selection involves a compromise. In this 
case, a Compromise must be made between computational cost, accuracy, and ease 
of implementation. In this chapter, I will discuss the various methods used in this 
study. The level of theory and basis set are discussed first, followed by the crystal- 
field simulation method. Finally, I will discuss the convergence criteria used in this 
study. 

LEVEL OF THEORY AND BASIS SET SELECTION 

The ultimate goal of quantum chemistry is the calculation of the solution to 
the Schrodinger equation resulting in an exact molecular wave function. However, 
that requires an infinitely large basis set expansion of the wave function along with 
the implementation of full configuration interaction. This is clearly an unattainable 
goal. In actuality, the basis set is truncated to a finite size and the electron 
correlation, if treated at all, is limited to a few configurations. These specifications 
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define a theoretical model within which all structures, energies, and other physical 
properties can be studied.!9 While proper selection of the model is important, 
recognizing and understanding the limitations of the given model on the results is 
imperative. 

The main sources of error in any geometry optimization using the gradient 
method are neglect of electron correlation ( where it is either ignored or 
incompletely treated) and the use of a finite basis set. The combined effect of these 
errors on the results obtained from a computed molecular structure are shown in 
Figure 2-1. The figure depicts the error in some arbitrary structural parameter 
(such as a C-H bond) in a variety of molecular environments plotted against 
increasing basis set size. When an infinite basis set is used (the Hartree-Fock 
limit), the error due to the neglect of electron correlation is found to be remarkably 
constant for the given parameter over a wide range of molecules. For much smaller 
basis sets, the error is expected to fall within the shaded area. However, at some 
point (X), the basis set is sufficiently large that the scatter in the calculated 
parameter 1s extremely small. This results in a constant known error for the given 
parameter called the offset value. The remaining absolute error is due to the 
truncation of the basis set and the neglect of electron correlation. Thus, the major 
flaw with small basis sets is the inconsistency of the error.* 

In order to keep the computational costs low for large molecules, the 
calculation should be performed as close as possible to the point X. Any increase 
in the basis set results in an increase in computation time and disk storage space 
proportional to the fourth power of the number of basis functions (N4). Beyond 


point X, the wave function approaches the Hartree-Fock limit where only electron 
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correlation error remains. The incorporation of any high-level electron correlation 
treatments results in a large increase in computation time and disk storage roughly 
proportional to n® where n is the number of electrons in the molecule.!! Even the 
most -caioait correlation methods, such as the local correlation method!2, are two 
to three times more expensive than a single determinant calculation. Therefore, 
depending upon the size of the molecule and the number of basis functions, the 
treatment of electron correlation may be prohibitively expensive for a given 
molecule. 

Based on the above observation, calculations can be done efficiently at point 
X. For many structural parameters, point X corresponds to approximately double 
zeta basis sets. However, for some parameters such as torsional angles around 
nitrogen or oxygen, it is essential to add at least one set of polarization functions to 
the basis to reach that point. The offset value for most bond angles is zero using 
this basis set. However, bond distances have non-zero offset values at this level 
and for some types of bonds the basis set must go beyond double zeta.!3 A 
common misconception is that the addition of d orbitals (polarization functions) to 
first and second row elements is incorrect since d orbitals in these elements are not 
actually occupied. The addition of polarization functions is only done to allow 
greater flexibility in the description of the bonds. Their addition does not imply that 
d orbitals are actually involved in the formation of these bonds. 

Because of the size of the H2SQO4 molecule, the treatment of electron 
correlation is prohibitively expensive. Unless stated otherwise, the calculations in 
my study were completed at the single determinant level using the gradient program 


TEXAS |4 with no symmetry restrictions. During the course of my study, Pulay's 





group completed a new version of TEXAS which was used for a few additional 
calculations. The major improvements are increased speed and the ability to use 
true Gaussian d functions. After some optimization of the program by Cordell for 
the Cray, I have seen up to a 30% increase in speed over the old version of 
TEXAS. The inclusion of true d functions has replaced the use of displaced p 
functions which were used to construct d functions in the previous version of 
TEXAS. Throughout this paper, the new version of TEXAS will be referred to as 
TEXAS(D). 

As stated earlier, one of the most important parts of any ab initio calculation 
is the selection of a basis set. The literature is inundated with almost 100 basis 
sets!5 which range from the minimal STO-3G basis set to very extensive triple zeta 
plus multiple polarization function basis sets. Thus, the selection of a basis set for 
a particular problem must be made based on the ab initio program used, the 
accuracy required, and how much computation time is available. The TEXAS 
program efficiently uses split valence shell basis sets where the s and p Gaussian 
exponents are equal. Pople's 6-31G** basis sets are used almost exclusively 
throughout this study.!6,17,18,19,20 Some calculations, which will be specifically 
mentioned, were completed using a 4-31G** basis set on sulfur2! with a 6-31G** 
basis set on oxygen and hydrogen (referred to throughout the remainder of this 
paper as 4-31G**). The asterisks in the notation for the basis sets refer to the 
inclusion of polarization functions. The first asterisk refers to the addition of d 
functions to the heavy atoms sulfur (exponent=0.65)!7 and oxygen 
(exponent=0.8)22. The second asterisk refers to the addition of p functions to 


hydrogen (exponent=1.1)22. As discussed earlier, the inclusion of polarization 
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functions is essential to accurately obtain torsional angles around oxygen. Another 
basis set which was used to study the sulfuric acid dimer was Pople's STO-3G*.23 
That basis set consists of the minimal STO-3G basis with the addition of five d 
functions only to second row atoms. 
METHOD TO SIMULATE CRYSTAL-FIELD 

Numerous studies have been made in order to apply ab initio calculations to 
solids. One reason is the advantages which ab initio calculations offer over some 
experimental methods. While X-ray and neutron diffraction provide valuable 
information about the spatial arrangement of the atoms in molecules and the packing 
of the molecules in lattices, they tell us nothing of the forces which determine the 
crystal structure. In addition, they give no information on why the molecules pack 
in the observed space group or the influence of crystal forces on the conformation 
of flexible molecules.24 In some cases the experimental methods cannot "locate" 
particular parts of a molecule, especially hydrogens. In the case of sulfuric acid, 
three X-ray diffraction studies, with the latest being done in 1978, could not 
"locate" the hydrogens. It was not until a neutron diffraction study was done in 
1983 that the hydrogens were located. Another problem with the experimental data 
is the vibrational averaging which is inherent to all experimental methods. It is 
these questions and others to which ab initio calculations can be directly applied. 

One area in which theoretical methods have been frequently applied is the 
study of hydrogen-bonding. The ab initio methods which have been used range 
from the complex use of an atom-centered multipole expansion29 to the simple 
point charge model.? It was not until the study by Saebg, et al. that the influence of 


crystal forces on the intramolecular geometry was studied in detail. Saebg reported 
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impressive agreement between X-ray diffraction and calculated structural 
parameters by using point charges based on a Mulliken population analysis. 
Regardless of the methods which have been used, the importance of the 
electrostatic interaction energy in determining the energetically most favorable 
geometry is well established.8.9:24.26.27 In fact, it led Bonaccorsi, et al.28 to the 
“electrostatic assumption". Bonaccorsi and coworkers postulated that the minima 
in the total interaction energy coincide with those in the electrostatic energy for 
rotational degrees of freedom.28 Thus, it is clear that approximating the 
electrostatic energy via a model can reveal valuable molecular structure information 
if the model is sufficiently detailed. 

One study which specifically considered crystal-field effects was the 
calculation of the rotational conformation of acetamide.2? In that case, the crystal 
structure had shown that the torsion angle was between 0° and 30°. However, 
several attempts by other authors to compute the lattice energies and total 
crystallographic conformational energies (lattice+torsional) could not account for 
that torsional angle. A procedure was developed to compute the lattice energy as 
the sum of three long-range contributions (electrostatic, polarization, and 
dispersion) and a short-range repulsive contribution. That procedure yielded a 
minimum lattice energy with a torsional angle between 0° and 30° as observed in the 
crystal structure.29 Similarly, through the use of calculated electrostatic molecular 
potentials, researchers were able to compute the preferred arrangement of linear 
chains of formamide molecules as observed in the crystal.39 

The combination of experimental and theoretical methods has been used in 


several studies. One study of hydrogen bonding, crystal packing, and the effect of 








crystal forces on the molecular conformation of amides and carboxylic acids 
expounded upon the "marriage" of X-ray diffraction and ab initio calculations.24 
The goal of that study was to obtain a general force field for organic and biological 
molecules. In that study, the lattice energy was calculated as a sum of the Lennard- 
Jones potential and the Coulombic interaction. The researchers were able to use 
their procedure to explain the crystal-packing mode of formic and acetic acid. Their 
study clearly showed that the electrostatic interactions were the overriding factor in 
determining the particular packing occurring.4 In another "marriage", 
experimental electron densities were combined with calculations in order to study 
the role of Coulomb forces in the crystal packing of amides. In this case, the 
researchers were able to use the Coulomb interaction energies in order to determine 
the possible and preferred molecular packing modes. In addition, they were able to 
obtain reliable estimates of the van der Waals atom-atom potential parameters using 
the calculated Coulomb energies.2/ 

In one investigation, the electrostatic energy in the formic acid crystal was 
studied.26 The geometry of formic acid was not optimized; instead the researchers 
used the experimental gas and crystal geometries to conduct their study. In their 
study, two approaches were pursued - multipole and point charge lattice sums. In 
the multipole method, the electrostatic energy is expanded in a senes to yield the 
various multipole moments. That particular study used multipole moments up to 
the sixth moment. The electrostatic energy then becomes a sum of the interaction 
energy between the polar multipole moments. In the point charge method, the 
charge distribution is represented by a set of point charges which were obtained 


from a Mulliken population analysis. The electrostatic energy simply reduces to the 
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sum of the interaction energy between the point charges. Several important 
conclusions were made from the results. First, the electrostatic portion of the lattice 
energy of formic acid was both large in magnitude and varied strongly with the 
orientation of the molecule. Several orientations, including the one observed 
experimentally, were found to minimize the electrostatic energy. The location of 
these minima were dictated by interactions with nearest-neighbor molecules. 
Second, the point charge method agreed qualitatively with the multipole series for 
the dependence of the electrostatic energy on the molecular orientation. Finally and 
most importantly, the electrostatic component of the lattice energy was found to act 
as the driving force for molecular distortion upon crystallization.2© 

The conclusion which can be drawn from all of the above studies is that 
theoretical calculations can give extremely valuable information about a wide variety 
of crystal parameters. However, it 1s now time to specifically address the crystal 
parameter which is of interest in my study - molecular structure. A straightforward 
approach which has been applied to the calculation of molecular structure is to 
compute it directly in the solid. This involves selecting a model of an extended 
"supermolecule" that is sufficiently large to represent the condition existing inside 
the crystal. That same method has been previously used to study such phenomena 
as the adsorption of a gas on a surface, where the solid is represented by a small 
number of atoms. One way of implementing this method is to consider all the 
molecules within a unit cell as the "supermolecule" and then impose boundary 
conditions to correct for effects at the cell edge. However, as the size of the 
molecular system increases, the accuracy of the computed results suffers. Asa 


result, the “supermolecule" method has not been conducive to highly accurate 
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results for molecules of appreciable size.4 A recent study using that method has 
been completed by Cordell and Boggs on furan molecules.3! The calculations were 
done on furan dimers and trimers located in the spatial orientation found in the 
crystal. Their study was able to show which neighbor molecules stabilized the 
crystal structure and to what extent. The interaction energies of the various subsets 
of the unit cell were shown to be essentially additive, with the packing forces 
expected to decrease the calculated crystal stabilization by 5 to 10%. In addition, 
they concluded that the short C=C bond length reported for the crystal was an 
artifact of the X-ray diffraction experiment.3! However, as stated earlier, the 
“supermolecule” technique is only useful for relatively small molecules. There is 
also an appreciable cost associated with it which increases as N4. Thus, large basis 
set and/or large molecule calculations using the "supermolecule" concept are very 
expensive. 

Based on the previous discussion, there are major problems which afflict 
many of the methods which have been used to simulate crystal-field effects. All of 
the above methods suffer from one or more of the following problems: high 
computational cost, low accuracy, or difficulty in implementation. Recently, a 
method has been proposed by Saebg, et al. to calculate crystal structures which 
seems to overcome these problems. Their study consisted of investigating the 
geometry changes which occur in cyanoformamide in going from the gaseous to the 
solid states, with emphasis on the effect of the crystal forces on the intramolecular 
geometry. Their approach consisted of using a simple point charge model to 
simulate the crystal-field effects. First, the geometry of the free molecule was 


optimized and a wave function obtained. Next, a Mulliken population analysis was 
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run to generate the magnitude of the point charges. The point charges were then 
placed around the free molecule in the spatial arrangement which corresponded to 
the previously determined X-ray diffraction structure. Since the hydrogens were 
not located in the X-ray diffraction experiment, the N-H bond distances were set to 
1.0 A. Finally, the geometry of the free molecule was optimized inside of the point 
charges. The resulting optimized structure yielded bond lengths and bond angles 
with a maximum deviation of 0.011 A and 0.4°, respectively, compared to the X- 
ray diffraction structure. The major advantages of their method were not only the 
accurate results and the ease of implementation but also the very slight increase in 
computation time compared to the free molecule. The small increase in computation 
time was only due to the increase in the number of one-electron integrals and 
nuclear repulsion terms.? Thus, their method appears to possess all the advantages 
of the previously discussed methods with none of the disadvantages. 

There are a few problems with the point charge model which need to be 
considered. The model assumes that the crystal-field effects can be simulated 
through purely electrostatic interactions. However, these are not the only forces 
which bind molecules together in solids. In order to fully understand the types of 
forces which bind crystals together, the various crystals which exist must be 
reviewed. There are basically four types of crystals - ionic, covalent, metallic, and 
molecular. Metals represent a special class of solids and will not be discussed 
further. In ionic crystals, the atoms are bound together by largely electrostatic 
forces and there is no such thing as an individual molecule. Covalent solids are 
actually one huge molecule which is composed of covalently bonded atoms. The 


atoms are held together by sharing electrons with their neighbors. On the other 





hand, molecular crystals actually contain a collection of distinguishable molecules 
which are bound together. The molecules are held together mostly by dipole 
interactions, van der Waals’ forces and hydrogen bonds. In reality, very few 
crystals fall strictly into one of these categories. Instead, many crystals are held 
together by more than one of these forces. Not only are there numerous examples 
of crystals composed of both ionic and covalent bonds but even molecular crystals 
may contain a small amount of ionic or covalent bonding.32 But, all of these forces 
share a common denominator - they all depend upon charged particle interactions. 
In summary, the forces which bind crystals together are very complex and in 
certain cases may not be fully described by only electrostatic interactions. 

The foremost problem with the point charge model used by Saeb¢ and 
coworkers is that it does not specifically allow for intermolecular electron transfer. 
The neglect of electron transfer in their method could certainly degrade the accuracy 
of calculated crystal structures in which covalent bonding is a major force. 
However, since few covalent crystals exist, their method should give satisfactory 
results for a wide variety of other crystals. In fact, the study by Saebg indicated 
that hydrogen bonding was dominated by electrostatic effects. While the neglect of 
electron transfer must not be forgotten when interpreting results using this method, 
the point charge method should produce reasonable results when applied to a large 
variety of crystals, especially molecular crystals involving hydrogen bonding. 

There still remain a couple of problems with the Mulliken population 
analysis which need to be considered. First, it tends to produce charges which are 
rather arbitrary and basis set dependent. Any defects in the point charge model are 


dominated by possible errors in the magnitude of the point charges.? Second, since 





18 


itis based on a simplified model of describing the electron distribution, the 
Mulliken population analysis often yields multipole moments which are quite 
different than those calculated from the actual wave function.® Although, there are 
other standard population analysis methods such as that of Léwdin, they suffer 
from the same disadvantages. ! 

In order to obtain more accurate charges, a different method has been 
recently proposed.® In that method, the molecular electrostatic potential evaluated 
at points in space around the molecule is used as a guide and is then fit to the point 
charge models. The procedure provides for selected restraints on the calculated 
dipole or quadrupole moments, and for the use of additional point charges to 
represent lone pairs in the molecule. That study found very good agreement with 
experimental enthalpies of formation for hydrogen bonding on nucleic acid pairs. 
Based on the above results, the researchers found their method superior to the 
Mulliken population analysis. Even given a set of point charges which accurately 
reproduce the molecular electrostatic potentials, it is not a tnvial task to determine a 
set of nonbonded parameters that can be used with the charges to give good 
geometries.8 This method also tends to assign larger absolute values of charge to 
the atoms compared to the Mulliken population analysis. However, ease of 
implementation was the overriding factor in choosing between these two methods. 
While their method may produce a more accurate point charge model than the 
Mulliken method, it is by no means a trivial task to implement their method. The 
method which provides the best compromise based on the selection criteria is the 
point charge model using the Mulliken population analysis. 


Both versions of TEXAS have been modified to accommodate the point 





19 


charge model. The programs require that the point charges be input in Cartesian 
coordinates with their associated net charge. The energy which is calculated 
excludes the point charge-point charge interactions but includes the point charge- 
molecule interactions. Therefore, the calculated energy accurately shows the 
stabilization of the crystal when compared to the free molecule energy. 

The major reason that the Mulliken population analysis is easy to implement 
is because the required parameters are already calculated in any standard SCF 
program. The only required parameters are the atomic number, the density matrix 
(P), and the overlap matrix (S). The number of electrons that are associated with 
any basis function is simply the diagonal element of the product of the density and 
overlap matrices which corresponds to that basis function. The corresponding 
number of electrons that are associated with a given atom in a molecule is obtained 
by summing over all basis functions on that atom. This assumes that all the basis 
functions are centered on atomic nuclei. The net charge of any atom (point charge 
magnitude) is simply the difference between the charge of the atomic nucleus and 
the number of electrons associated with that atom. The net charge is positive if the 
number of electrons associated with the atom is less than the nuclear charge and 
negative otherwise. Although this method is simple, it must be remembered that 
there is no unique definition of the number of electrons that are associated with a 
given atom or nucleus in a molecule. In the case of CHy4, the Mulliken population 
analysis results in a net charge on hydrogen which is double that given by the 
Léwdin method at the STO-3G level ( 0.06 and 0.03, respectively ). The basis set 
used also has a profound effect on the net charge which results. A Mulliken 


population analysis on CHy4 gives hydrogen a net charge of 0.06 at the STO-3G 
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level but 0.12 at the 6-31G** level.! Despite its drawbacks, the Mulliken 
population analysis is readily obtained and has been shown to give exceptional 
results. 

As discussed previously, the method employed by Saebg, et al. was simply 
to optimize the free molecule surrounded by the point charges. They chose to place 
the point charges in positions as determined by X-ray diffraction and to assign them 
net charges based on a Mulliken population analysis of the optimized free molecule. 
Despite their convincing results, they suggested two improvements within the 
frame of their model. The first was to introduce the new net charges calculated for 
the free molecule in the surrounding point charges. Secondly, they recommended 
increasing the size of the basis set.? These questions and several others are the 
focus of my study. 

Based upon the suggestions of Saeb¢ and coworkers, the expansion of the 
basis set by inclusion of polarization functions is a logical extension of their work. 
The increased flexibility offered by the polarization functions is an improvement in 
the basis set over that used in their study. Also, the idea of introducing the new 
point charge values calculated from the free molecule suggests an iterative approach 
to me. Although Saebg and coworkers envisioned only one additional optimization 
cycle, I felt that multiple optimization cycles would allow more flexibility in the 
model by allowing it to dynamically change. The iterative technique was applied in 
my study by optimizing the free molecule inside the point charges based on both X- 
ray and neutron diffraction structures. Then the net charges were recalculated by 
running a Mulliken population analysis on the new "free" molecule geometry. The 


new net charges were assigned to the experimental point charge locations and the 
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cycle was repeated until the free molecule converged in both structure and energy. 
This should eliminate some of the error associated with the Mulliken population 
analysis. 

One problem with the Saeb¢ method is that it requires a detailed 
experimental structure of the solid. While X-ray diffraction structures are available 
for a wide range of compounds, they quite often are unable to locate the hydrogens 
which was the case for both cyanoformamide and sulfuric acid. Then, in order to 
implement their method, the positions of the hydrogens must be assumed which 
certainly could have profound effects on the results. In addition, experimental 
results are prone to mistakes just like any other method. For sulfuric acid, the early 
X-ray diffraction results determined the symmetry to be Cay instead of the actual 
C2. Therefore, it would be helpful to have a method which was not quite so 
dependent on the experimental results. 

One solution is to position the point charges based upon the free molecule 
instead of the experimental structures. The coordinates of the free molecule can 
easily be converted from Cartesian to fractional coordinates and then used to 
generate the point charges. That requires only the cell dimensions and the space 
group of the crystal and it preserves the internal geometry of the free molecule. 
These parameters can be determined from experiment, assumed from similar 
crystals, or simply guessed. This would permit calculations on molecules for 
which there are no experimental results or for which incomplete crystal structures 
exist. However, the major advantage of my method would be for calculations on 
molecules where the electron density is shifted away from the nuclear center and for 


which only X-ray diffraction crystal structures exist. One of the disadvantages of 
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X-ray crystallography is the fact that it gives bond lengths which are not very 
accurate. That 1s due to the manner in which X-rays interact with matter. The 
interatomic distances determined by X-ray diffraction represent the distance 
between centers of electron density and not nuclei. In molecules of this type, X-ray 
diffraction can give results for bond lengths which greatly differ from the true bond 
lengths.’ The error in bond length could cause profound effects on the structure 
calculated in the point charge model. The technique of basing the point charges on 
the free molecule has been applied in my study using both the X-ray and neutron 
diffraction data. In these calculations, the only experimental information that was 
used were the cell dimensions and the space group. The calculations were repeated 
in an iterative manner until the "free" molecule converged in both structure and 
energy. By basing the point charges on the "free" molecules, any errors caused by 
inaccurate X-ray bond lengths should be eliminated. 

Since the value of net charge determined by the Mulliken population 
analysis seems to be somewhat arbitrary, the structural results should have a 
dependence on the charges that are generated. However, it is not known whether 
this is a significant factor affecting the calculated structure. The dependence on the 
magnitude of the charges was approached in two manners. First, a Mulliken 
population analysis was run on both the experimental X-ray and neutron diffraction 
structures. That required a single point SCF calculation be done at the experimental 
geometry in order to generate a wave function. The resulting point charges varied 
by as much as 15% from those calculated for the optimized free molecule. Second, 
the value of the point charges was normalized to one of the previously calculated 


charges. This preserved the neutrality of the molecule and resulted in charges for 





2 


the atoms about double those for the free molecule. The new charges were than 
used to compute geometries which were compared to those calculated using the 
charges based on the free molecule. 

Although the increased computation time required by the addition of point 
charges is small, on the order of 500-700 seconds per run, it would still be 
beneficial to limit the point charges to as few as possible. The increased 
computation time becomes a significant factor for the iterative methods which 
require up to 30 runs. In the Saebg study, the closest neighbor molecules were 
included. But there are obviously pieces of these molecules which are far enough 
away from the free molecule such that their electrostatic interaction is essentially 
zero. Hence, if they were not included, the computation time could be further 
reduced. In order to investigate the effect on the computed structure, calculations 
were made where only atoms within a given sphere from the center were included. 
Because of the requirement to maintain a neutral charge, a dummy charge was 
generated to make the sum of the point charges zero and it was placed in excess of 
150 A from the free molecule. 

CONVERGENCE CRITERIA 

The point charge calculations require a different type of convergence criteria 
than normal calculations on isolated molecules. This is due to the residual 
Cartesian forces which result from the presence of the point charges. In the point 
charge calculations, the optimized geometry occurs when the internal forces are 
small but the Cartesian forces may still be quite large. 

An isolated molecule geometry is considered optimized when the following 


conditions are met:33 
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. All Cartesian forces are less than 0.005 mdyne/angstrom. 


Stretching and bending internal coordinate forces are less than 
0.005 mdyne/angstrom or /radian. 

Torsional internal coordinate forces are less than 0.001 
mdyne/radian. 


The internal force on torsional coordinates has changed sign. 


Conditions 1 through 3 are the normal conditions used to determine 


geometry convergence. Condition 4 is absolutely essential to ensure convergence 


of low energy torsional and out-of-plane motions.24 The above conditions should 


give convergence of calculated bond lengths to + 0.002 A, bond angles to + 0.2°, 


and torsional angles to + 2.0°. 


As discussed earlier, the convergence criteria for the point charge 


calculations must be slightly different. For these calculations, conditions 2 through 


4 above must be met with the additional requirement for all Cartesian coordinate 


displacements to be less than 0.0005 A. The restriction on the change in Cartesian 


coordinates is based upon the work by Saebg and my review of other isolated 


molecule optimizations when conditions 1 through 4 were met. 





CHAPTER THREE 
SULFURIC ACID 

One reason for the great interest in sulfuric acid is because of its presence in 
“acid” rain. The harmful effects of "acid" rain have been widely reported over the 
last decade. These detrimental effects include large fish kills, contamination of 
water supplies, and damage to automobile finishes. The Northeastern portion of 
the United States has been hit especially hard. Research has found this to be 
predominantly due to the increased burning of high sulfur coal in the Midwest. 
More recently, the effects of "acid" rain have been reported worldwide. Many old 
churches and monuments in Europe are rapidly being destroyed by "acid" rain. For 
these reasons and several others, such as the detection of sulfuric acid in the 
Venusian atmosphere, many studies have been done on sulfuric acid. These 
studies have centered on the formation and structure of sulfuric acid, especially the 
intermediates that are involved. Results from some of these studies have lead to an 
understanding of the formation of sulfuric acid and its detection in the vapor phase. 

As discussed in Chapter One, sulfuric acid undergoes significant structural 
changes transforming from a gas to a solid. Regardless of the phase, the symmetry 
remains C7. The H2SO4 molecule with C2 symmetry is shown in Figure 3-1 along 
with the atom designations that will be used throughout the paper. These structural 
alterations can easily be seen in a comparison of the gas and solid structures as 
shown in Table 3-1. The most dramatic change is the shortening of the S-O single 
bond on the order of 0.05 A. In addition, there is a noticeable decrease in the 
O=S=O bond angle coupled with an increase in the O-S-O bond angle. These 
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Figure 3-1 
Atom Numbering Scheme for Sulfuric Acid (H2SQOg4) 
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Table 3-14 
Sulfuric Acid (H2SO4) Gas and Crystal Geometries 


MwWb = X-RAY(1978)° ND4 X-RAY(1965)® 


Phase gas solid solid solid 
S-O1 1.574+0.01 1.528(S5) 1.48 1.535(15) 
S=02 1.422+0.01 1.419(5) 1.49 1.426(15) 
O1-H O77 5001 if 0.91 f 
ZH-O1-S NG Sot 1) f 113.6 f 
ZO1-S-O1’ lpi saall 103.7(4) 112.8 104.0(10) 
Z02-S-O2' Way Saal 118.3(4) 110.1 118.6(10) 
ZO1-S-O2 108.6+0.5 110.6(3) 108.0 110.5(10) 
ZO1-S-O2' 106.4+0.5 106.3(3) 109.0 105.9(10) 
t(H-O1-S-O1')8 90-1 f 84.1 ie 
t(H-O1-S-O02)8 -20.8+1 fi -36.5 ie 


a. Bond lengths in angstroms. Bond angles in degrees. 

b. Preferred microwave structural parameters and uncertainties from reference 35. 

c. X-ray diffraction geometry with estimated standard deviations (esd's) in 
parentheses from reference 36. 

d. Neutron-diffraction 10°K structure from reference 37. The estimated 
uncertainties are reported in the paper as +0.05-0.10 A. 

e. X-ray diffraction geometry with esd's in parentheses from reference 38. 

f. The hydrogens were not located in either X-ray diffraction study. 


g. T(ijkl) is the torsional angle between the ijk and jkI planes. 
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significant changes are attributed to the hydrogen bond formation in the crystal.’ I 
have been able to reproduce these structural changes with the point charge model. 

The remainder of this chapter will deal with the results I have obtained for 
both the gas and solid phases of sulfuric acid. The gas phase will be discussed 
first, followed by the crystal phase. 

GAS 

Although vapor pressure studies in 1923 showed that solutions with 
concentrations exceeding 85% had measurable vapor pressures above 200° F, early 
studies of gaseous sulfuric acid were hampered by several problems. An infrared 
(IR) spectrum of gaseous sulfuric acid was not completed until 1965, because of its 
low volatility and attack on windows at high temperatures. In that study, the 
observed H2SO4 bands were assigned in agreement with the expected chemical 
series (X-SO2-Y: X,Y= F, OH, Cl, .. .,CH3). Based on those assignments, the 
researchers incorrectly concluded that free H2SO4 molecules had an approximate 
tetrahedral configuration with Coy symmetry similar to the other X-SOQ2-Y 
compounds.3? 

The first molecular orbital calculation of the gas phase structure of H2SO4 
was undertaken in 1978. The purpose of that study was to gain insight into the gas 
phase reaction of SO3 and H20 to form H2SQz4 using semi-empirical CNDO/2 
calculations. Those researchers incorrectly calculated the symmetry of H2SO4 as 
Cv in agreement with the earlier IR study. The calculated structural parameters 
were O1-H= 1.04 A, S-O1= 1.63 A, S-O2= 1.56 A, ZO1-S-O1'= 97°, ZH-O1-S= 
97.5°, and ZO2-S-O2'= 116° with the H-O1-S-O1'-H' atoms being planar.40 


Their configuration is referred to as the compact Coy geometry where the OH 
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groups point toward and bisect the terminal SO? group with a torsional angle (ZH- 
O1-S-O1' and ZH'-O1'-S-O1) of 180°. 

The first ab initio calculation on H2SO4 was completed in 1980 due to the 
importance of sulfur-oxygen compounds in air pollution chemistry, specifically 
"acid" rain.4! The optimum geometries in that study were calculated using the 
STO-3G* basis set. Previous calculations on similar compounds using the STO- 
3G or 4-31G basis sets had shown poor agreement between experiment and 
calculated structural parameters though the STO-3G* basis set had been shown to 
give acceptable geometries. Since the researcher's principal interest was in the 
geometry about the sulfur atom, the OH length and SOH angle were fixed at 0.994 
A and 105°, respectively [these values corresponded to their optimum values in 
SO2(OH)]. The calculated STO-3G* geometry for sulfuric acid was S-O1= 1.63 
A, S-02=1.45 A, Z01-S-O1'= 99°, ZO1-S-O2= 107°, and 202-S-O2'= 125°. 
The H-O-S-O dihedral angles were set to 90° by the researchers based upon the best 
idealized conformation. The addition of d functions to the basis set was shown to 
greatly improve the calculated energy. The energies calculated at the optimum 
geometry for the STO-3G* and 4-31G* basis sets showed an energy improvement 
of 316.1 and 170.6 kcal/mole, respectively, over the energies at the STO-3G and 4- 
31G levels.4! 

The correct C2 symmetry was finally determined in 1980 from the 
microwave spectra of four isotopic species of H2SO4. From those spectra, the 
researchers determined a detailed structure for gaseous H2SO4 which is shown in 
Table 3-1. In contrast to the previously determined compact Coy geometry, the OH 


bonds were found to have rotated past the SO2 and SO2' bonds resulting in 
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dihedral angles of -20.8° with them. That structure possessed a H-O1-S-O1' 
dihedral angle of 90.9° as compared to 180° found in the CNDO/2 calculation.35 
The H2SQ04 molecule with Co symmetry is shown in Figure 3-1. 

The continued interest in H2SO4 prompted two additional ab initio studies 
of its structure and torsional modes.42,43 Both of these studies determined the 
potential energy surface for H2SQO4 as a function of 6 (ZH-O1-S-OL') and 9' (Z 
H'-Ol'-S-O1). The p,p structure ("p" denotes periplanar or cis) is shown in 
Figure 3-2 where 6=6'=0°. These researchers performed most of their calculations 
at the STO-3G level except for some selected conformations which were calculated 
at the 4-31G level. Despite the different structures used for sulfuric acid, as 
discussed below, the results of these two studies were surprisingly similar. 

The first study was started prior to the completion of the microwave study. 
As a result, Kaliannan and coworkers used the geometry of dimethyl] sulfate 
[SO2(OCHS3)] obtained by electron diffraction in their study of sulfuric acid. The 
O-H bond length and H-O-S bond angle were set at 0.96 A and 109.47°, 
respectively. Their construction of the potential energy surface consisted of 
calculations at 22 conformations of } and $'. They found the global minimum 
occurring at d='=90° in good agreement with the microwave study value of 
o=o'=90.9°. Their calculated dipole moment at that conformation was 2.51 debye 
(D) compared to the microwave value of 2.725 D. A secondary minimum was 
found to occur at d=-6'=120° with an energy about 1.2 kcal/mole above the global 
minimum. They also determined that the variation in energy was very small in the 
range 90°-120° for both and 9’, suggesting that the molecule is reasonably flexible 


in that range. A comparison of their potential energy surface to that of H2POq- 





Figure 3-2 


H2SQgq p,p structure 6=0'=0° ("p" denotes periplanar or cis) 
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found the two to be very similar with the global minimum occurring at 0=0'=90° in 
both cases. The results of their 4-31G calculations at the STO-3G optimized 
geometries confirmed the presence of the global minimum at 6=d'=90° in agreement 
with the microwave experiment.42 

Lohr has also conducted a study of the structure and torsional modes of 
H2SQj4. In his study, the microwave values for bond lengths and bond angles 
were used. He constructed his potential surface from 43 calculated points at the 
STO-3G level. Lohr found the global minimum had C2 symmetry and occurred at 
d=6'=98.6°, somewhat larger than the microwave value of d=o'=90.9° . A local 
minimum with Cs; symmetry was also found at d=-d'=104° with an energy of 1.36 
kcal/mole relative to the global minimum. Calculations at eleven selected points 
were done using the larger 4-31G basis set. The structure with C2 symmetry was 
found to have a minimum energy at d='=94.9° which is somewhat closer to the 
microwave value than is the STO-3G value.43 The results of both torsional mode 
studies confirm the C2 symmetry which is observed in the microwave study. 

I optimized the gas phase of sulfuric acid with several different basis sets. 
The programs TEXAS, TEXAS(D), and GAUSSIAN-82 “4 were used to obtain an 
optimum structure at the 6-31G** level. In addition, calculations were conducted 
at the 4-31G** level using TEXAS and at the STO-3G* level using GAUSSIAN- 
85 45. The results of my calculations on gaseous sulfuric acid are presented in 
Table 3-2. 

The 6-31G** basis set admirably reproduces the experimental structure. 
With very few exceptions, the calculated 6-31G** structure is within the 


experimental uncertainties of the microwave study. The major variances occur in 





S-Ol 


S=02 


O1-H 


ZH-O1-S 


Z201-S-O1' 


ZO02-S-O2' 


Z01-S-O2 


ZO1-S-O2' 


t(H-O1-S-O1') 


t(H-O1-S-O2) 


Energy: 


Hartree 


kcal/mole 


a. Bond lengths in angstroms. Bond angles in degrees. 


Table 3-24 


Sulfuric Acid (H2S0O4) Gas Phase Geometries 


Mwb 


1.574+0.01 
1.422+0.01 
0.97+0.01 
108.5+1.5 
101.3+1 
13,31 
108.6+0.5 
106.4+0.5 
ot | 


-20.84+1 


6-31G**¢ 


LO 


1.411 


0.950 


110.7 


101.8 


1235 


108.1 


106.7 


88.9 


-23.2 


-698.05 1422 


-438,033.781 


6-31G**¢ 

TEXAS(D) 
1.569 
1.411 
0.950 
110.8 
101.8 
123.6 
108.1 
106.7 
89.2 


-22.9 


-698.052805 


-438,034.649 


6-31G**°¢ 


G-82 


1.569 


1.411 


0.950 


110.8 


101.8 


123.6 


108.1 


106.6 


88.8 


=23:3 


-698.052804 


-438,034.65 


4-31G** 


1.563 


1.408 


0.950 


LEZ 


101.7 


1235 


108.3 


106.5 


87.8 


-24.2 


-697.659157 


-437,787.631 
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SiO-3G7° 
G-85 
1.620 
1.446 
i993 
105.6 
100.1 
125.6 
108.9 
105.2 

Ten 


-31.0 


-689.766048 


-432,834.63 


b. Preferred microwave structural parameters and uncertainties from reference 35. 


c. This work. Calculations run with TEXAS unless noted otherwise. 
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the torsional angles which are notoriously difficult to reproduce. Considering the + 
2.0° uncertainty in the value of the calculated torsion angles, the agreement is very 
good. 

The development of TEXAS(D) by Pulay's group prompted me to optimize 
the gas phase structure with it for comparison. As shown in Table 3-2, the 6- 
31G** structures are identical for both TEXAS and TEXAS(D). My calculations 
confirm that the displaced p functions used to construct d functions in TEXAS are 
more than adequate. The only difference between these two calculated structures is 
in the energy. The TEXAS(D) energy is 0.868 kcal/mole lower than that calculated 
in TEXAS. The improved energy in TEXAS(D) must be due to the use of true d 
functions. For most calculations, this energy difference is insignificant. Based on 
calculated structures and energy, neither TEXAS nor TEXAS(D) offer any 
advantage over the other with the exception of the slightly lower energy in 
TEXAS(D). In terms of computation time, TEXAS(D) is superior, posting 
decreases of up to 30% over TEXAS. TEXAS(D), however, requires double the 
amount of disk space to store the integrals as compared to TEXAS. 

Because of the recent acquisition of a Cray version of GAUSSIAN-82, I 
decided to make a comparison between it and TEXAS(D). Since both programs 
incorporate true d functions, I expected identical results. As shown in Table 3-2, 
the results are essentially identical. The few tenths of a degree difference between 
torsional angles is well within the calculated uncertainty of + 2.0°. The TEXAS(D) 
program calculated an energy which is only 6 x 10-4 kcal/mole lower than 
GAUSSIAN-82. Based upon my comparison, both TEXAS(D) and GAUSSIAN- 


82 do give identical results for calculations conducted with the same basis set. 





Initially, I was interested in reducing the computation time required for one 
optimization run. As a result, I completed an optimization run of gaseous sulfuric 
acid at the 4-31G** level. The 4-31G** basis set I utilized contained a 4-31G** 
basis on sulfur and 6-31G** on both oxygen and hydrogen. As shown in Table 3- 
2, the 4-31G** basis set gave structural parameters almost identical to those 
computed with TEXAS at the 6-31G** level. The major distinctions are the shorter 
O1-S bond length and the torsional angles at the 4-31G** level. The computation 
time saved at the 4-31G** level was significant. The 4-31G** calculation ran 15% 
faster than the 6-31G** calculation, 1761 seconds compared to 2066 seconds. The 
time savings is due to the reduction in the number of primitive Gaussian functions. 
Although the number of contracted Gaussians is the same for both basis sets at 89, 
the number of primitive Gaussians is reduced from 208 to 198 1n going from the 6- 
31G** basis set to the 4-31G** basis set. The major disagreement between the 
two calculations is the total energy. The 6-31G** calculation results in an energy 
247 kcal/mole lower than the energy at the 4-31G** level. For many calculations, 
the 4-31G** basis set would be satisfactory. Since I was concerned with 
observing structural changes, I chose to use the more accurate 6-3 1G*™* basis set in 
the crystal calculations. In addition, the 6-31G** basis set {(11s4p1d/4s 1p) 
contracted to [4s2p1d/2s1p]}} represents a significant improvement over the basis 
set used by Saebg and coworkers {(7s3p/4s) [4s2p/2s]}. The 6-31G** basis set 
not only increases the number of primitive functions but also adds polarization 
functions to hydrogen and the heavy atoms. 

In order to make comparisons with some dimer calculations, I optimized 


sulfuric acid at the STO-3G* level. As shown in Table 3-2, the calculated 
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parameters at this level differ greatly from the experimental structure. The large 
error in the torsional angles at the STO-3G* level confirms the need for polarization 
functions on oxygen to accurately calculate these angles. More importantly, while 
no symmetry restrictions were placed on any of the other calculations, I had to 
restrict the STO-3G* calculation to C) in order to prevent the program from 
determining a structure with Cy; symmetry. While minimal basis set calculations 
can yield valuable information, the need for accurate structural information almost 
always requires using a larger basis set. 
CRYSTAL 

The first X-ray diffraction study of sulfuric acid was completed in 1954 by 
Pascard. That study incorrectly concluded that sulfuric acid crystallized in the 
noncentrosymmetric space group C2 and possessed Coy symmetry.*© In 1964, 
Pascard-Billy refined the original intensity data with anisotropic temperature factors 
for oxygen atoms to a reliability (R) factor of only 0.124. The reported dimensions 
of the H2SO4 molecule deviated significantly from idealized Coy symmetry. 
Pascard-Billy determined that the correct space group was actually C2/c and the 
H2S04 molecule possessed C2 symmetry. Although Pascard was not able to locate 
the hydrogens, he theorized that they each formed a hydrogen bond with a double 
bonded oxygen on an adjacent molecule. In that configuration, one H2SO4 
molecule is hydrogen bonded to four adjacent molecules with the hydrogen atom 
lying on the -O1...02= line.38 

In 1978, Yu and Mak felt that the atomic parameters reported earlier could 
be further refined by full-matrix anisotropic least squares with correction for 


secondary extinction. Their study showed significant improvement of the standard 
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deviations in atomic parameters over those reported earlier with an R value of 
0.099. Their measured hydrogen bond length (O1-H...O2) of 2.624 A belongs to 
the short type. They concluded that the hydrogen atom almost lay on the 
internuclear O1...02 line if a tetrahedral S-O-H angle was assumed.3® 

The hydrogen atoms were finally located in 1983 by Moodenbaugh and 
coworkers using powder neutron-diffraction techniques. Diffraction data was 
collected near the melting point (240°K) and at 10°K to search for order-disorder 
transitions in the hydrogen bond network. The neutron-diffraction scans revealed, 
to the limits of detection, that the sulfuric acid sample was single phase with 
monoclinic symmetry at both temperatures. The principal conclusions from their 
study were that the hydrogens were ordered at both temperatures and the 
parameters of the other atoms were qualitatively consistent with the single-crystal 
X-ray studies.3/ 

The generation of the point charges used in my calculations was facilitated 
by the use of several computer programs. The BMFIT*7 (Best Molecular Fit) 
program derives the best least-squares fit between two sets of atoms. I used 
BMFIT in order to convert calculated structures into fractional coordinates in the 
crystal space. The program DAESD‘8 calculates interatomic distances and angles 
in crystal structures. The DAESD program was used to generate all the molecules 
within a given unit cell. The programs ZERO and MOLE were written by me in 
order to generate the actual point charge locations. Both programs convert the 
contents of a unit cell into a molecule and its nearest neighbors. In the program 
ZERO, all neighbor atoms within a given sphere of radius r from the central 


molecule are punched. In addition, a dummy charge is generated and placed at 





(100,100,100) in order to maintain the neutrality of the point charges. In MOLE, 
all neighbor molecules that have their central atom within r from the central 
molecule are punched. In the case of MOLE, the neutrality of charge is maintained 
by including whole molecules instead of individual atoms. 

My studies have shown that great care must be taken in generating the point 
charges. As discussed in Chapter Two, I initially tried to limit the number of point 
charges in order to reduce the computation time. I did this by using the program 
ZERO to generate the point charges. The program ZERO allowed me to include 
only the nearest atoms and not the nearest molecules. I could then eliminate atoms 
on adjacent molecules that were far enough away that they had essentially zero 
electrostatic interactions with the central molecule. However, this method proved 
to be entirely unsatisfactory. Not only was the savings in computation time 
negligible, but this method also introduced artificial symmetries into the calculated 
structure. In my particular case, the point charges generated using the program 
ZERO artificially reduced the symmetry of H2SOq4 from C2 to C;. As a result, all 
of my crystal calculations were made using the nearest neighbor molecules 
generated by MOLE. After careful scrutiny of the point charges, I found that 
including the 18 nearest H2SO4 molecules insured that all nearest neighbor atoms 
were included in the point charges. All of the reported calculated crystal structures 
used 126 point charges based on the 18 nearest H2SO4 molecules. 

The crystal calculations which I have completed can be divided into four 
categories - SF, SE, MF, and IF. These designations are used throughout the 


remainder of the paper and are defined as follows: 
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A single optimization cycle was run. (Note: An optimization cycle 
means that the convergence criteria of Chapter Two have been met.) 
The point charge magnitudes are based on the optimized free (gas 
phase) molecule. The point charge locations are based on the 
experimental crystal structures. This is the method used by Saeb¢ 
and coworkers. 

This is identical to the SF method except that the point charge 
magnitudes are based on a calculation at the experimental crystal 
geometry. 

This method is a continuation of the SF method involving multiple 
optimization cycles. In this case, a Mulliken population analysis is 
run on the new "free" molecule geometry obtained from the SF 
method. (Note: The term "free", in quotation marks, refers to the 
central molecule surrounded by the point charges.) The new net 
charges are assigned to the experimental point charge locations and 
another optimization cycle is completed on the "free" molecule. The 
process is repeated until the "free" molecule converges in both 
structure and energy. 

This is an iterative method involving multiple optimization cycles 
similar to the MF method. In this case, however, both the point 
charge magnitudes and locations are based on the optimized free 
(gas phase) molecule. The only information used from the 
experimental crystal study is the unit cell dimensions and space 


group. After each optimization cycle, the magnitude and location of 





the point charges are adjusted to those corresponding to the newly 
optimized "free" molecule. The process is continued until the "free" 
molecule converges in both structure and energy. 

The SE method was conceived by me to test the sensitivity of the calculated 
crystal structures to the magnitude of the point charges. A Mulliken population 
analysis was run on the wave function from a single point SCF calculation at the 
experimental geometry to obtain the experimental net charges. As expected, the net 
charges from the experimental geometry differed significantly from those of the 
optimized free molecule. Although using the net charges from the free molecule 
seems logical since a wave function already exists, no study has been made 
previously to compare the effects of the two sets of net charges. Since a single 
point SCF calculation is relatively small (about 650 seconds), the calculation of the 
experimental net charges and their subsequent use in the SE method seems a 
reasonable alternative to the SF method. 

The use of the iterative method in the point charge model is an extension of 
an improvement proposed by Saeb¢g and coworkers. They simply suggested that 
the point charge model might be improved by including the net charges of the 
newly optimized "free" molecule. My idea was to iterate this process until the 
molecule converged in both structure and energy. My iterative technique is 
therefore self-consistent. The major advantage of my method 1s that it adds 
flexibility to the point charge model. A significant concern with the point charge 
model is that it may not incorporate all the forces which bind crystals together. As 
discussed in Chapter Two, all of the binding forces in crystals involve charged 


particle interactions. Since the model basically allows for electrostatic interactions 
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between the electrons and the point charges, the iterative method should accurately 
recreate the environment inside the crystal. In addition, any effect of the arbitrary 
assignment of net charge in the population analysis should be greatly minimized 
since the net charges are varied in response to their environment. 

I devised the IF method for a couple of reasons. First, experimental crystal 
data is often incomplete or inaccurate. As discussed previously, X-ray 
crystallography frequently cannot locate hydrogen atoms. Even when all the atoms 
are located, the precision with which an atom can be positioned in a crystal is 
subject to systematic errors. These errors can be divided into two classes: errors in 
the model and errors in the data. If it was possible to collect data with no errors, 
many significant structural errors could be introduced if an improper model is used 
in the refinement. The IF method is not as dependent on the experimental data 
since it uses only the cell dimensions and space group of the crystal. The cell 
dimensions and space group are more easily determined and can be found even if 
all the atoms are not located. Second, a principal concern with the point charge 
model is that it may not adequately represent the crystal environment. The IF 
method allows for movement of the point charges. Therefore, the point charges 
should ultimately reside in the crystal locations if the point charges properly model 
the crystal environment. Thus, the IF method serves as a test of the validity of the 
point charge model. 

The first calculations that I completed were based on the X-ray study of Yu 
and Mak. Since the hydrogens were not located in their study, I assigned their 
location by assuming an O1-H bond length of 0.95 A and a hydrogen bond angle 
(ZO1-H...02)of 180°. The results of these calculations are presented in Table 3-3. 
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In addition to the calculations based on the X-ray study, I also conducted 
calculations based on the neutron-diffraction study. The 10°K neutron-diffraction 
data was used as the basis for the calculations because of its minimal thermal 
motion and due to errors in the 240°K data I discovered. The results of these 
calculations are presented in Table 3-4. The numerous problems with the neutron- 
diffraction study will be discussed later. 

One aspect that all of the calculated solid structures share is a lower energy 
than the gas phase. The energy improvement in the solid structures is on the order 
of 50-70 kcal/mole. The lower energy in the crystal is expected based on the 
stabilizing effect of the neighbor molecules in the solid. The calculated energy 
improvement which corresponds to 12.5-17.5 kcal/mole per hydrogen bond 
compares favorably to 7.8 kcal/mole for cyanoformamide?, 5.4 and 5.9 kcal/mole 
for CH3CN...HOCH3, respectively4?, 3.9-9.6 kcal/mole for the water dimer 
(depending upon basis set)°9, and 7.0-10.3 kcal/mole for the various 
monohydrated formamide complexes?!. In addition, all of the calculated crystal 
structures mimic the major structural changes which occur in H2SO«g in going from 
a gas to a solid. Specifically, the calculated crystal structures show a substantial 
shortening of the S-O bond along with an increase in the ZO1-S-Ol' and a 
decrease in the ZO2-S-O2' compared to the gas. Further, all of the calculated 
crystal structures agree with the general and very characteristic structural feature 
between bond angles for sulphones (XSO2Y). Hargittai has demonstrated that the 
ZX-S-Y < ZX(Y)-S=O < Z O=S=O, regardless of the X and Y ligands. The 
Z O=S=O is considerably larger than the ideal tetrahedral angle while the Z X-S-Y 


is considerably smaller.52 While these results alone substantiate the point charge 
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Table 3-34 
Sulfuric Acid (H2SO4) X-ray Crystal Geometry 
XRAY(1978)> SFC SEC MFe IF¢ 
S-O1 1.528(5) 1.558 oot 1eS53 1552 
S=O02 1.419(5) 1.422 1.424 1.428 1.428 
O1-H 0.9504 0.971 0.972 0.978 0.978 
O1...02¢ 2.624(8) 2.553 DSBS 2.550 D556 
ZH-O1-S 109.94 113.9 iba 114.7 114.8 
Z01-S-O1' 103.7(4) 104.3 104.6 105.3 LOS 
20228202 118.3(4) 120.8 120.5 119.6 119.4 
Z@OiES-O2 110.6(3) 107.3 107.3 Oe 107.4 
Z01-S-O2' 106.3(3) 108.0 108.1 108.4 108.3 
ZS-O1...02®  110.0(3) =—111.3 111.3 111.3 UII 
ZO1...02°-S®  157.9(4) 160.7 160.7 160.5 160.4 
t(H-O1-S-Ol') — 86.34 92.3 92.0 91.6 89.5 
t(H-O1-S-O2) ——-27.34 Oe, -22.6 B2ou/ -25.7 
Energy: 
Hartree -698.130620 -698.155672 
-698.13468 1 -698.161332 
kcal/mole -438,083.479 -438,099.199 
-438,086.027 -438, 102.750 


a. Bond lengths in angstroms. Bond angles in degrees. 


b. X-ray diffraction geometry with esd's in parentheses from reference 36. 
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This work. Calculated using the 6-31G** basis set. 

. Based on setting the O1-H bond length to 0.95 A and the O1-H...02 bond 
angle to 180°. 

Atom is located at fractional coordinates -1/2+x,1/2+y,z on an adjacent 


molecule. 


ts 





45 


Table 3-44 
Sulfuric Acid (H2SO4) Neutron-Diffraction 10°K Geometry 
ND> SF¢ SE¢ MF IFC 
S-O1 1.48 1.554 e551 1.549 esiay2 
S=02 1.49 1.425 1.428 1.431 1.428 
O1-H 0.91 OES 0.976 0.980 0.979 
O1...024 acai 255 ZS Zo) 255 
H...024 1.65 NS, 1.58 1.58 1.58 
ZH-O1-S 136 HSS 114.4 114.3 NS 
ZQO1-S-O1' 112.8 104.9 105.7 106.0 105.4 
Z02-S-O2' LORI 120.0 119.1 118.6 eS 
ZO1-S-O2 108.0 Uy Oy 2 107.0 1073 
ZO1-S-O2' 109.0 108.3 108.4 108.8 NOSE 
ZO1-H...024 lS Wy Nes) 17358 Ls i227 
ZS-O1...024 P93 7211 112.0 112.1 Wa) 
ZO1...029-S4 155.3 oi),8 159.1 Ihsrs)| See 
t(H-O1-S-O1') 84.1 oom 90.5 92-6 91.3 
t(H-O1-S-O2) -36.5 -21.9 -25.1 -23.4 -23.9 
Energy: 
Hartree -698.136173 -698.165534 
-698.160937 -698.160931 
kcal/mole -438 086.963 -438, 105.387 
-438,102.503 -438, 102.499 
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Bond lengths in angstroms. Bond angles in degrees. 


. Neutron diffraction 10°K structure from reference 37. The estimated 


uncertainties are reported in the paper as +0.05-0.10 A. 
This work. Calculated using the 6-31G** basis set. 
. Atom is located at fractional coordinates -1/2+x,1/2+y,z on an adjacent 


molecule. 
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model, the remaining differences between the calculated and experimental crystal 
structures need closer scrutiny. 

The differences between the calculated and neutron-diffraction crystal 
structures are easily explained by examining the experimental data. Moodenbaugh 
and coworkers report estimated standard deviations (esd's) of + 0.05-0.10 A. 
Although their calculated esd's are about 0.01 A, the larger esd's arise from an 
apparent discrepancy in the O1-H bond length which is known to be temperature 
independent to about + 0.002 A.37 Based on their + 0.05-0.10 A uncertainty in 
bond lengths, all of the calculated values fall within that range. Any comparison 
with the neutron-diffraction study by Moodenbaugh and coworkers is essentially 
useless because of the extremely large uncertainties in the experimental data. 

There are numerous reasons for the large errors reported by Moodenbaugh 
and coworkers. First, the structural model that is used in any refinement of the raw 
crystal data is crucial to the results. In their study, they based the heavy atom 
positions on the 1965 X-ray study with its larger esd's instead of the 1978 study 
with its smaller esd's. They were unaware of the 1978 study. They also chose to 
uSe isotropic temperature factors for simplicity, even though both of the X-ray 
diffraction studies indicate that anisotropic factors would be more appropriate. 
Their R factors, on the order of 0.20, are approximately double those from the X- 
ray diffraction studies. Second, the quality of their data is rather poor with an 
unfavorable typical peak to background ratio of 1:2. This is probably due to the 
fact that they used a powder sample since they could not obtain a suitable single 
crystal. The high background caused by incoherent scattering required them to use 


the Rietveld refinement method which is not normally used. 
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Finally, the data from the neutron-diffraction experiment is questionable at 
best. Table 3-5 shows the bond lengths and angles reported by them at 240°K. 
However, when I used their coordinates to calculate structural parameters, the 
actual S-O1 bond length differs greatly from that they reported. At first glance, 
only the reported S-O1 bond length would appear to be in error. If the value of 
1.67 A that I calculated is correct, the S-O1 bond length would undergo a 
shortening of 0.19 A between 240°K and 10°K which is highly unlikely since 
H2S 0a was found to be ordered at both temperatures. In addition, the values I 
calculate from their data for the ZO1-S-Ol' (=87.9°) and ZO02-S-O2' (=133.1°) 
differ greatly from those in the X-ray diffraction study (103.7° and 118.3°) and in 
the neutron-diffraction study at 10°K (112.8° and 110.1°). The discrepancies in the 
S-Q1 bond length and the angles about sulfur suggest to me that there is an error in 
the fractional coordinates of the sulfur atom. Since sulfur occupies the special 
position O,y, 1/4, the error appears to be in the y fractional coordinate. I have made 
several attempts to determine the correct sulfur coordinates by varying the y 
coordinate. In my opinion, the most probable structure for H2SOq, based on the 
data presented in the paper at 240°K, occurs when y=0.068 as shown in Table 3-5. 
The reported value is y=0.028(6). My value of 0.068 appears to be a logical 
transposition of the reported value. At y=0.068, there is much better agreement 
between the 10°K and 240°K structures. I have talked to Moodenbaugh about this 
problem. Moodenbaugh was not previously aware of the problem and told me he 
put much more reliance on the heavy atom positions from the X-ray studies. 
Further, he indicated absolutely no interest in trying to resolve the problem. The 


numerous problems with the neutron-diffraction study by Moodenbaugh and 
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Table 3-54 
H2S0O4 Neutron-Diffraction Structural Parameters at 240°K 
Reported> Calculated Probable 

S-O1 1.54 1.67 1.54 
S02 1.40 1.40 1.49 
O1-H 1.00 1.00 1.00 
Oi 2.47 2.47 2.47 
..J@2¢ 1.49 1.50 Ie 
ZO1-H...02¢ 163 164 164 
Z0O1-S-O1' 87.9 97.9 
Z02-S-O2' Neieis 119.6 


a. Bond lengths in angstroms. Bond angles in degrees. 

b. Neutron diffraction 10°K structure from reference 37. The estimated 
uncertainties are reported in the paper as +0.05-0.10 A. 

c. Atom is located at fractional coordinates -1/2+x,1/2+y,z on an adjacent 


molecule. 
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coworkers certainly warrant another neutron-diffraction study of H2SQq. 

A comparison of the calculated and X-ray crystal structures in Table 3-3 
shows relatively good agreement. By far the largest discrepancy is in the S-O1 
single bond length. The X-ray experiment determined the S-O1 distance to be 
about 0.03 A shorter than the calculation. As discussed earlier, a major 
disadvantage of X-ray diffraction is that it determines distances between centers of 
electron density and not nuclei. Since the shorter bond length could be a 
consequence of the X-ray experiment, I decided to investigate it further. In order to 
investigate the S-O1 bond length discrepancy, I chose three distinct approaches. 
The first was to calculate the crystal structure with an expanded basis set. The 
second involved calculations on the H2SQO4 dimer. In the third method, I used the 
program MOCALC:S in order to look at the total electron density along the S-O1 
bond. 

Since the point charge model may not adequately describe the hydrogen 
bonding, the calculated S-O1 bond length may be too long. However, various 
studies on dimers have shown that not only similar double zeta basis sets4991,54,55 
but also even the minimal STO-3G basis set°° adequately describe hydrogen 
bonding. Saebg and coworkers concluded from their study of cyanoformamide 
that hydrogen bonding was dominated by electrostatic effects.? In addition, Tang 
and Fu concluded from ab initio studies of hydrogen bonding between methyl 
cyanide (CH3CN) or methyl isocyanide (CH3NC) and methanol (CH30H) that it 
was a result of ionic binding of two closed-shell systems.4? Nevertheless, in order 
to see if an expansion in the basis set might better describe the hydrogen bonding, a 


calculation was made with the 6-31 +G** basis set.!0 The "+" implies the addition 
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of a diffuse sp shell to the hydrogen atoms (exponent=0.0360). The addition of 
diffuse functions is required in the calculation of anions, transition structures, and 
other species with significant electron density far removed from the nuclear center. 
The optimized geometry with the 6-31 +G** basis set was practically identical to 
that of the 6-31G** basis set with no difference in the S-O1 bond length. 

Another approach I pursued was to optimize the structure of the H2SO4 
dimer. The two molecules chosen were those that were hydrogen bonded in the 
crystal. The dimer calculation showed a definite shortening of the S-O1 bond at the 
STO-3G* level. The hydrogen bonded S-O1 and S-O2 bond lengths were 1.58 A 
and 1.46 A, respectively, compared to 1.62 A and 1.45 A in the monomer at the 
STO-3G* level. The shortening shown in the dimer must be viewed with caution. 
In the optimized dimer, the two molecules rotated, as a result of the program, such 
that two hydrogen bonds were formed between them instead of only one as found 
in the crystal. This is apparently caused by the lack of neighbor molecules which 
prevent the rotation in the crystal. My attempts to conduct dimer calculations at the 
6-31G** level with the inclusion of the point charges have not been successful. I 
initially planned to optimize the dimer structure at the STO-3G* level with the 
inclusion of the point charges. However, because of the success of the total 
electron density plot in explaining the shortened bond in the X-ray study, I have not 
completed any additional dimer calculations. 

The program MOCALC?? allows plotting electron density of either a single 
molecular orbital or all of the molecular orbitals. MOCALC was specifically 
written to use the data generated in the TEXAS program. In order to investigate the 


shortened S-O1 bond, I have completed a total electron density plot of H2SO4 in 





the O1-S-O1' plane which is shown in Figure 3-3. As Figure 3-3 clearly shows, 
the electron density of the Ol atom is definitely shifted toward the S atom. After 
talking with several crystallographers, I have determined that there is no commonly 
used quantitative method to account for this apparent bond shortening. The 
crystallographers simply accept the fact that it exists. In order to get a quantitative 
feel for the magnitude of this shortening, I have simply assumed the location of the 
center of the electron density on the O1 atom in Figure 3-3. By making simple 
measurements with a ruler, I have determined, that in the case of H2SQq, the 
shortening of the S-O1 bond length in the X-ray study is on the order of 0.03-0.05 
A. Thus, as shown in Figure 3-3, I conclude that the short S-O1 bond length 
reported in the X-ray studies 1s a consequence of the X-ray experiment. 

The various methods used in this study to determine the magnitude and 
location of the point charges yielded some interesting results. As shown in Table 
3-3 and 3-4, the SF and SE methods yielded almost identical geometries despite the 
differences in the magnitude of the point charges. In both cases, the magnitude of 
the point charges obtained from the experimental crystal structures in the SE 
method differed by as much as 15% from those calculated for the optimized free 
molecule. The differences between the SF and SE structures were no more than 
0.004 A in bond lengths and 0.9° in bond angles. In most cases, within the 
uncertainties of the calculations, the two methods gave structures which agree. The 
convincing agreement between these two methods shows that the structures 
calculated using the point charge method are relatively insensitive to the magnitude 
of the charges. This indicates that the arbitrary manner in which the Mulliken 


population analysis assigns net charge has little effect on the calculated results. 





Figure 3-3 
Total Electron Density Plot of H2SOq in the O1-S-O1' Plane 
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My study shows that the iterative technique used in the MF and IF methods 
offers several advantages over the single optimization cycle methods. Both of these 
methods show substantial improvements in the total energy over the SF method. 
The energy improvements are on the order of 15-19 kcal/mole. Throughout the MF 
and IF optimization cycles, the calculated structural parameters continued to better 
approximate those in the crystal. Specifically, the S-O1 bond length continued to 
shorten, the ZO1-S-Ol1' increased, and the ZO2-S-O2' decreased. The 
improvements in both the energy and structural parameters confirm the use of the 
iterative approach. One particularly interesting result from the MF and IF methods 
is the change in gross atomic populations. Three ab initio studies by Ottersen and 
Jensen found that there is a transfer of electron density from the hydrogen atom to 
the acceptor atom. The formamide (HCONH?2) dimer showed a transfer of 0.05 
electron across the hydrogen bond.*4 Ottersen's study of the keto (HCONH2) and 
enol (HOCHNH) tautomers of formamide exhibited gross atomic population 
increases of 0.057 and 0.089 on the acceptor atoms and decreases of 0.046 and 
0.044, respectively, in the bridging hydrogens.°> Their study of the four 
monohydrated formamide complexes demonstrated gross atomic population 
increases of 0.009-0.034 for the acceptor atoms and decreases of 0.015-0.039 for 
the bridging hydrogens.>! The gross atomic populations for my calculations using 
the iterative method are presented in Table 3-6. All four of my calculations clearly 
follow the trend of Ottersen and Jensen by exhibiting a transfer of electron density 
from the bridging hydrogens to the acceptor atoms (O2/02'). These results 
confirm the ability of the point charge model to accurately describe the crystal 


environment for hydrogen bonded systems. In summary, the iterative technique 
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Table 3-64 
H2S04 Gross Atomic Populations 
Free Optimized Net 
Molecule Geometry? Change 
X-ray: 
MF- S 14.18843 14.11860 -0.070 
O1/OL 8.65419 8.67702 +0.023 
H/H’ 0.60258 0.51039 -0.092 
02/02’ 8.64902 8.75329 +0.104 
IF- S 14.18843 14.10879 -0.080 
O1/O1 8.65419 8.67743 +0.023 
H/H'’ 0.60258 0.50813 -0.094 
02/02’ 8.64902 8.76004 +0.111 
Neutron-Diffraction: 
MF- S 14.18843 14.11909 -0.069 
O1/Ol’ 8.65419 8.66920 +0.015 
H/H’ 0.60258 0.50504 -0.098 
02/02’ 8.64902 8.76621 +0.117 
[F- S 14.18843 14.10930 -0.079 
O1/O1' 8.65419 8.67745 +0.023 
/Et 0.60258 0.50761 -0.095 
02/02' 8.64902 8.76029 +0.111 
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a. All values were calculated from the Mulliken population analysis with the 6- 
31G** basis set. 

b. These values were obtained using the indicated methods after the calculations 
converged in both energy and structure. In all cases this required 5 


optimization cycles. 
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tends to minimize any detrimental effects of the Mulliken population analysis. In 
addition, by allowing the point charge model to dynamically change in relation to its 
environment, any effects or forces which are not entirely electrostatic by nature 
appear to be accounted for. The only disadvantage of the iterative technique is the 
amount of computation time required. With each run averaging 2,500 seconds and 
up to 30 runs required, the iterative technique is fairly expensive. However, the 
accurate results from these methods far outweigh their cost. The MF and IF 
methods are clearly superior to the SF method used by Saebg and coworkers. 

My study has shown that the IF method 1s particularly useful. Not only 
does it allow calculations on molecules for which incomplete experimental 
structures exist, but it also gave better results in my calculations based on the X-ray 
data. The fact that the MF and IF X-ray structures are almost identical substantiates 
the point charge model. The point charge model reproduces the deformation forces 
which are present in the solid upon crystallization as shown by the IF method. In 
the IF method, the point charges are initially positioned in the free molecule 
geometry but they are subsequently adjusted to the newly calculated "free" molecule 
geometry. Since the point charges eventually assume the approximate structure of 
the crystal, the point charge model must essentially reproduce the forces present in 
the crystal. Based on my results, both the MF and IF methods are superior to the 
SF method and they allow much more flexibility in the model. Since both the MF 
and IF methods have been shown to give almost identical results, either method can 
be used but the IF method should certainly be used in cases where the experimental 
data is incomplete. 


A major concern with the point charge model is the dependence of the 
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calculated structure on the magnitude of the point charges. In order to investigate 
the dependance, I made several runs in which I varied only the magnitude of the 
point charges. The point charges were positioned according to the arrangement 
found in the 10°K neutron-diffraction study and the free molecule was optimized 
inside these point charges as in the SF method. The results of those calculations 
are shown in Table 3-7. In run #1 which uses the SE method, the magnitude of the 
point charges was calculated from the experimental geometry. The net charges in 
run #2 were calculated from the optimized free molecule. Run #2 is equivalent to 
the SF results reported in Table 3-4. In run #3, I arbitrarily normalized the net 
charges of the free molecule to the O2/O02' net charge. Despite the difference in the 
net charges which vary by up to 80%, the calculated geometries are extremely 
similar. The maximum difference in bond lengths and bond angles is 0.02 A and 
2.1°, respectively, occurring only in run #3 with its extreme values of charge. I 
conclude from these results that the structures calculated using the point charge 
model are fairly insensitive to the magnitude of the point charges. Therefore, the 
arbitrary assignment of net charge by the Mulliken population analysis has little 


effect of the calculated structure. 





Table 3-74 
Calculated H2SO4 Solid Phase Geometries 


#16 #26 #36 
Net Charge: 
S 1.75069 Ie Sills y 2.79124 
O1/OL’ -0.56293 -0.65419 -1.00796 
H/H 032572 0.39742 0.61234 
O2/02' -0.70833 -0.64902 - 1.00000 
Structure: 
S-O1 1551 1.554 1.549 
S=02 1.428 1.425 1.436 
O1-H 0.976 0.972 Or 2 
ZH-O1-S 114.4 LS ys: [S26 
ZQO1-S-OL' Gt 104.9 106.6 
ZO02-S-O2' ed 120.0 118.1 
ZO1-S-O2 O72 107-2 106.4 
ZO1-S-O2' 108.4 1O8.3 109.4 
t(H-O1-S-O1') 90.5 O31 oles 
t(H-O1-S-O2) -25.1 -21.9 -25.6 
Energy: 
Hartree -698.160937 -698.136173 -698.194785 
kcal/mole -438,102.503 -438,086.963 -438,123.743 
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Bond lengths in angstroms. Bond angles in degrees. 
. Geometry optimized at the 6-31G** level. The point charge locations are based 
on the 10°K neutron-diffraction data. The magnitude of the point charges were 


obtained as explained in the text. 





Chapter Four 
CONCLUSION 

In this paper, I have presented a thorough ab initio study of sulfuric acid 
and the point charge model. The point charge model has been shown to be 
extremely advantageous in the calculation of crystal structures with few limitations. 
In addition, I have presented a comprehensive review of various crystal field 
simulation methods which have been used. In this chapter, the important results 
from my study are summarized and my suggestions for future work are offered. 

My calculations of the gas phase structure of H2SO4 show that the 6-31G** 
basis set accurately reproduces the microwave structure. With few exceptions, the 
calculated 6-31G** structure 1s within the experimental uncertainties of the 
microwave study. The 4-31G** basis set gave a similar structure to the 6-31G** 
results with a substantial savings in computation time. My study has confirmed 
that the use of displaced p functions to simulate d functions has essentially no effect 
on the calculated structure. The only difference being in the total energy which is 
less than 1 kcal/mole lower with the true d functions. The ab initio programs 
TEXAS(D) and Gaussian-82 give identical results with the 6-31G** basis set for 
H2S0q4 as expected. My calculations confirm that polarization functions are 
absolutely essential to reproduce torsional angles on oxygen. 

My calculations of the crystal structure of H2SO4 demonstrate that the point 
charge model accurately reproduces the structural trends which occur in 
transforming from the gas to solid phase. My study of the point charge model 
proves that it precisely simulates the deformation forces which are present in the 


60 





61 


solid upon crystallization. During the course of my study, I discovered that the 
only neutron-diffraction study of H2SOz is plagued by numerous errors, rendering 
it totally useless. 

The calculated crystal structures based on the X-ray diffraction study are in 
good agreement with the experiment. The major discrepancy is in the S-O single 
bond length which is calculated about 0.03 A longer than in the X-ray study. After 
a thorough investigation, I conclude that the discrepancy is due to a shifting of 
electron density from the oxygen atom towards the sulfur atom, resulting in the 
short experimental bond length. 

The iterative technique gives superior results to any of the single 
optimization cycle methods. Even with the point charges arranged in the gas 
geometry, the point charge model gives almost identical structural results to those 
calculations where the charges are positioned in the crystal structure. The point 
charge model more than adequately describes the crystal forces, especially for 
systems involving hydrogen bonding. 

The iterative technique minimizes any possible drawbacks of the point 
charge model. Any consequence of the arbitrary assignment of net charge by the 
Mulliken population analysis is minimized by allowing the charges to vary. By 
allowing the point charge model to dynamically change in relation to its 
environment, any forces which are not electrostatic by nature appear to be 
accounted for. 

Further evidence of the point charge model's ability to represent the crystal 
environment is exhibited in the change in gross atomic population throughout the 


iterative cycles. There is a definite shift of electron density from the bridging 
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hydrogens to the acceptor atoms (O2/02'), identical to those found in other ab initio 
studies of dimers. 

My study proves that the calculated crystal structures are insensitive to the 
magnitude of the point charges. I conclude that the arbitrary assignment of charge 
by the Mulliken population analysis has a negligible effect on the point charge 
model. 

SUGGESTIONS FOR FUTURE WORK 

The remainder of this chapter is a short presentation of areas that I feel need 
further work. 

Although the results of my study confirm the validity of the point charge 
model, the only types of crystals that have been tested to date involve hydrogen 
bonding. I strongly feel that further studies should be done on crystals bound 
together by both ionic and covalent bonds. 

A quantitative method needs to be developed to treat the shift of electron 
density that frequently occurs in X-ray diffraction studies. Currently, no 
quantitative method exists and bond lengths determined by X-ray diffraction are 
often reported shorter than their actual value. The development of such a method is 
desperately needed to facilitate comparisons between calculated and experimental 
crystal structures. 

The problems that I have presented with the neutron-diffraction study of 
H2S0q4 by Moodenbaugh and coworkers point to repeating that study. The new 
study should incorporate the Yu and Mak data and be done on a single crystal if 
possible. As a minimum, Moodenbaugh’'s data should be refined again using the 


heavy atom positions from the 1978 X-ray study and anisotropic temperature 
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factors. I expect that the results of a new neutron-diffraction study will confirm 
that the S-O single bond length is closer to my calculated value than the X-ray 
value. 

After discussions with several crystallographers, they all feel that a present 
day X-ray diffraction study of H2SO4 could locate the hydrogens. With the aid of 
a crystallographer, I have unsuccessfully tried to refine the Yu and Mak data to 
locate the hydrogens. Since no further knowledge can be gained from past X-ray 
diffraction studies, I feel a new X-ray study should be completed. The positions of 
the hydrogens, if located, can then be compared to the neutron-diffraction study to 


establish its accuracy. 
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